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Abstract 

In this paper, for the first time, a three-dimensional treatment of microtubules' polymerization is pre- 
sented. Starting from fundamental biochemical reactions during microtubule's assembly and disassembly 
processes, we systematically derive a nonlinear system of equations that determines the dynamics of mi- 
crotubules in three dimensions. We found that the dynamics of a microtubule is mathematically expressed 
via a cubic-quintic nonlinear Schrodinger (NLS) equation. We show that in 3D a vortex filament, a generic 
solution of the NLS equation, exhibits linear growth/shrinkage in time as well as temporal fluctuations 
about some mean value which is qualitatively similar to the dynamic instability of microtubules. By solving 
equations numerically, we have found spatio-temporal patterns consistent with experimental observations. 



* Corresponding author: Division of Experimental Oncology, Cross Cancer Institute, 11560 University Avenue, Ed- 
monton, AB T6G 1Z2 Canada, Tel: +1(780)432 8906, Fax: +1(780)432 8892; Email: jtus@phys.ualberta.ca 



1 



I. INTRODUCTION 



More than 25 years ago, Del Giudice et al. [1] argued that a quantum field theory approach 
to the collective behaviors of biological systems is not only applicable but the most adequate as it 
leads naturally to nonlinear, emergent behavior which is characteristic of biological organization. 
They followed a line of reasoning championed by Davydov [2, 3] and Frohlich [4] who emphasized 
integration of both conservative and dissipative mechanisms in biological matter leading to the emer- 
gence of spatio-temporal coherence with various specific manifestations such as almost lossless energy 
transport and long-range coordination. 

Microtubules (MTs) are long protein polymers present in almost all living cells. They participate 
in a host of cellular functions with great specificity and spatio-temporal organization. In most multi- 
cellular organisms, the interior of each cell is spanned by a dynamic network of molecular fibers called 
the cytoskeleton ('skeleton of the cell'). The cytoskeleton gives a cell its shape, acts as a conveyor 
for molecular transport, and organizes the segregation of chromosomes during cell division, amongst 
many other activities. The complexity and specificity of its functions has given rise to the notion that 
along with its structural and mechanical roles, the cytoskeleton also acts as an information processor 
[5], or simply put the "cell's nervous system" [6]. A microtubule is a hollow cylinder, a rolled-up 
hexagonal array of tubulin dimers arranged in chains along the cylinder ('protofilaments'). Within 
cells, microtubules come in bundles held together by 'microtubule associated proteins' (MAPs). The 
geometry, behavior and exact constitution of microtubules varies between cells and between species, 
but an especially stable form of a microtubule runs down the interior of the axons of human neurons. 
Conventional neuroscience at present ascribes no computational role to them, but models exist in 
which they interact with the membrane's action potential [7, 8]. 

Microtubules are very dynamic bio-polymers that simply lengthen and/or shorten repeatedly at 
the macroscopic level during a course of time on the scale of minutes. At the microscopic level, 
however, several biochemical reactions are taking place in order for an individual microtubule to 
undergo an assembly or disassembly process. This dynamical behavior of microtubules (so-called 
dynamic instability) has attracted many investigators for decades to examine microtubules' behavior 
in many aspects (see section 2 for further details). 

Although there is no systematic description for the microtubule's assembly /disassembly process 
at the microscopic level, several theoretical models have been proposed to describe the macroscopic 
or statistical aspect of lengthening/ shortening of microtubules using nonlinear classical equations [9- 
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14]. The latter studies while providing good agreements with experimental results, they are more or 
less phenomenological. Therefore, several features of the microtubule's assembly /disassembly process 
might not be adequately captured. 

In this paper, we propose a systematic model for the microtubule's assembly/disassembly process 
at the microscopic level using a first-principles quantum mechanical approach as a starting point. 
In this model we consider an individual microtubule with length L consisting of iV tubulin layers 
viewed here as a quantum state \N). The state can be raised/lowered by a creation/annihilation 
operator (i.e. polymerization/depolymerization process) to the \N + 1)/\N — 1) state. The corre- 
sponding microtubule is then longer /shorter by one tubulin layer from the original one. Based on 
the chemical binding reactions that are taking place during microtubule polymerization, a quan- 
tum mechanical Hamiltonian for the system is proposed. Equations of motion are then derived and 
transformed from the purely quantum mechanical description to a semi-classical picture using the 
inverse Fourier transformation. The resulting nonlinear field dynamics reduces to the cubic-quintic 
nonlinear Schrodinger equation that provides a richer dynamics than the previous phenomenologi- 
cal descriptions and includes both localized energy transfer and oscillatory solutions both of which 
have either been experimentally demonstrated or theoretically predicted earlier. We believe that 
this treatment is both useful and necessary to address the fundamental issues about the observed 
dynamical behavior of MTs. As stated by Del Giudice et al. [1] : "Systems with collective modes 
are naturally described by field theories. Furthermore, quantum theory has proven to be the only 
successful tool for describing atoms, molecules and their interactions." 

II. MICROTUBULE ASSEMBLY BACKGROUND 

A very rigid (by biological standards) and typically several micrometers long rod-like polymer plays 
an essential role during cell division. The so-called microtubule (MT) is assembled by tubulin poly- 
merization in a helical lattice. These protein polymers are responsible for several fundamental cellular 
processes, such as locomotion, morphogenesis, and reproduction [15]. It is also suggested that MTs 
are responsible for transferring energy across the cell, with little or no dissipation. 

Both in vivo and in vitro observations confirmed that an individual MT switches stochastically be- 
tween assembling and disassembling states that makes MTs highly dynamic structures [16, 17]. This 
behavior of MTs is referred to as dynamic instability. Dynamic instability of MTs is a nonequilib- 
rium process that has been the subject of extensive research for the past two decades. It is generally 
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believed that the instability starts from the hydrolysis of guanosine triphosphate (GTP) tubulin that 
follows by converting GTP to guanosine diphosphate (GDP). This reaction is exothermic and releases 
~ 8k B T energy per reaction [18], i.e. approximately 0.22 eV per molecule [19]. Here k B is the Boltz- 
mann's constant and T is the absolute temperature. Since GDP-bound tubulin favors dissociation, 
a MT enters the depolymerization phase as the advancing hydrolysis reaches the growing end of a 
MT. This phase transition is called a catastrophe. As a result of this transition, MTs start breaking 
down, releasing the GDP-tubulin into the solution. In the solution, however, reverse hydrolysis takes 
place and a polymerization phase of MTs begins utilizing assembly-competent tubulin dimers. The 
latter phase transition which comes after a catastrophe is called a rescue. Therefore, MTs constantly 
fluctuate between growth and shrinkage phases. 

Odd et al. [20] studied both experimentally and theoretically MT's assembly to extract their 
catastrophe kinetics. These authors proposed that a growing MT may remember its past phase 
states by analyzing growth characteristic of both plus and minus ends of several individual MTs. 
Their results showed that while the minus end growth time follows an exponential distribution, the 
plus end fits a gamma distribution. The exponential (gamma) distribution suggests a first (non-first) 
order transition between growing and shrinking phases. Statistically, the exponential distribution 
represents that the new state happens independently of the previous state. As a result, a MT with 
first order catastrophe kinetics does not remember for how long it has been growing. In contrast, 
the catastrophe frequency of a MT with non-first order kinetics would depend on its growth phase 
period. The gamma distribution suggests that the catastrophe frequency is close to zero at early 
times, increases over time and reaches asymptotically a plateau. This is consistent with observations 
that the catastrophe events are more likely at longer times. Odd et al. [20] concluded that such 
behavior implies that a 'crude form of memory' may be built in MT's dynamic instability. As a result, 
a microtubule would go through an 'intermediate state' before a catastrophe event takes place. 

The dynamics of transitions between growing and shrinking states is still a subject of controversy. 
It is suggested that a growing MT has a stabilizing cap of GTP tubulin at the end which keeps it 
from disassembling [16, 17]. Whenever MT loses its cap, it will undergo the shrinking state. Several 
theoretical and experimental studies have been devoted to the cap model. For the purpose of this 
paper, we emphasize the link between GTP hydrolysis and the switching process from growing to 
shrinking of a MT. GTP hydrolysis is a subtle biochemical process that carries a quantum of biological 
energy and thus allows us to make a link between quantum mechanics and polymer dynamics. We 
return to this theme later in the paper but first discuss the statistical methods used in this area. 
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A. Ensemble dynamics of microtubules 



As we discussed earlier, the MT dynamical instability has been the subject of numerous studies. 
Although the dynamical instability of MTs is a nonlinear and stochastic process, so far only their 
averaged behaviors have been analyzed using a simple model. Introducing p g (x, t) and p s (x, t) as the 
probability density of a growing and shrinking tip, respectively, of a MT with length x at time t, 
Dogterom and Leibler [14] proposed the following equations for the time evolution of an individual 
MT: 

dtPg = -fgsPg + fsgPs ~ V g d x p g , (1) 
dtPs = fgsPg ~ fsgPs ~ V s 8 x p s . (2) 

Here f gs and f sg are the transition rates from a growing to a shrinking state and vice versa. The 
average speeds of the MT in the assembly and disassembly states are given by v g and v s , respectively 
(see also [10, 11, 13]). 

Random fluctuations about the MT's tip location can be also modeled by adding a diffusive term 
in the above equations: 

dtPg = - fgsPg + fsgPs ~ V g 8 x p g + D g 8 xx p g , (3) 
dtPs = fgsPg - fsgPs - V s 8 x p s + D s d xx p s , (4) 

where D g and D s are the effective diffusion constants in the two states [21, 22]. 

Equations (3) and (4) describe the overall dynamics of an individual MT without considering the 
dynamics of GDP and GTP tubulin present in the solution. It is clear that the MTs are growing 
faster in the area with a higher concentration of GTP tubulin. Using this fact, Dogterom et al. [12] 
generalized the above model by incorporating the tubulin dynamics. They added two more equations 
to the above system: 

8 t c T = -v g s p g + kc D + L>iV 2 c T , (5) 
d t c D = v s s Ps - kc D + D 2 V 2 c D , (6) 

where ct and c D are average concentrations of GTP and GDP tubulin, respectively. D\ and D 2 are 
the diffusion coefficients, k is the rate constant and < So < 1. In view of the link to quantum 
transitions between GDP and GTP at the root of this problem we now introduce a method that allows 
a smooth transition from quantum to classical (nonlinear) dynamics of MT assembly /disassembly 
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process. Recently, Antal et al. [23, 24] also proposed a one- dimensional statistical model to describe 
the dynamic instability of MTs. They considered a two state model that MT grows with a rate A 
and shrinks with a rate ji and then obtained similar fluctuations in MT's length by varying A and \x 
as observed for MTs. 

The above mentioned studies generally considered a MT as a one-dimensional mathematical object 
that grows and shrinks with random rates. Such a treatment is not only applicable to MTs but also 
to a polymer system whose assembly or disassembly processes have an element of randomness. As a 
result, the biophysical and biochemical characteristics of MTs cannot be adequately captured. In this 
paper, however, we based our model on fundamental biochemical reactions that are occurring during 
microtubule's assembly and disassembly processes. This allows us to derive a nonlinear system of 
equations that determines the structure, the dynamics and the motion of MTs in 3D. 

B. Deriving Semiclassical Equations 

The underlying method we use here has been developed in a number of papers and a book [25-31] 
and is essentially semiclassical in nature. The treatment is quantitative in that important terms 
which are retained are calculated exactly and those which are very small but nevertheless significant 
are discussed at a later stage and their effect is estimated. The motivation for the method and 
a derivation of the dynamical field equation are presented in [27] and a discussion of the types of 
classical field solutions is presented in [28] . A fuller version has been published in the review paper in 
[31] whereas a very brief overview is given in [29]. It has been successfully applied to the phenomenon 
of superconductivity [26, 30] and when combined with topological arguments yields, for example, the 
correct temperature dependence of the critical current density in low temperature superconductors. 
One can also obtain the position of phase boundaries in metamagnets where previously only elaborate 
numerical techniques could provide this information [25]. Spatial correlations are fully incorporated 
using a renormalization technique and quantum fluctuations have been included also. It has been 
demonstrated that even when the method is generalized to include spin-dependent fields, the equation 
of motion for the field is of the same form [32] and the classical field equation is also of the same 
form for both Boson and Fermion particles. This does not mean that the Fermionic character of 
the electrons disappears because the statistics of the particles reappears in the choice of the classical 
field which satisfies the physical boundary conditions on the charge density. The method is basically 
non-relativistic although it could be readily generalized but here we use the non-relativistic version. 
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The starting point in this method is to write a generic form of second- quantized Hamiltonian using 
one particle state annihilation and creation operators: 

H = ^kSkSk + kA K 1, m9k9l9mgk+l-m, (7) 

k k, 1, m 

where the vectors k, 1 and m are shorthand labels for quantum numbers of a complete orthonormal 
set of particle functions in the usual way and we use the linear momentum conserving form for 
the two-body interaction. Depending on the system studied, using Fermi-Dirac or Bose-Einstein 
statistics one can derive the Heisenberg's equation of motion: 

ihd t q k (r,t) = [H,q k (r,t)\. (8) 

Now both sides of Eq. (8) are multiplied by Or 1 / 2 exp(— irj ■ r)a r) (t) and summed over r\. At the 
same time the matrix elements c^k and Ak, i, m are each expanded to second order in the deviations 
from the point (k , lo, m ). After a considerable amount of algebra and a series of transformations 
we find 

+ higher order terms, (9) 

where 

i/j(r,t) =fi- 1/2 ^exp(- ? r;-r)a T? (t), (10) 
v 

where Q is the volume over which the members of the plane wave basis are normalized [27]. Here ^ 
or /ij are constant parameters, determined by matrix elements and A kj \ m and their derivatives 
calculated at point (k , lo, m ). To convert Eq. (9) to a PDE in a complex number (c-number) field, 
rather than an operator, the center of expansion (k , lo, m ) is selected to be a critical or fixed point 
of the system. The reason for this is that close to a critical point it is an excellent approximation to 
replace the full quantum field, ip(r,t), by a classical component, ip c [33-35]: 

i>(r,t) =^ c (r,t)I + 0(r,t), (11) 

where I is the unit operator in Fock space, ip c is a c-number field, is a quantum mechanical operator 
with magnitude about |0| ~ h\ip c \ [36] (see [37] for details). 

In the next section we apply the above method to study the dynamic instability of an individual 
microtubule. 
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III. A QUANTUM MECHANICAL PICTURE OF THE MICROTUBULE ASSEMBLY 
PROCESSES 

A. Particle states 

Consider an individual microtubule in a free tubulin solution containing a large number of GTP- 
tubulin, GDP-tubulin and a pool of free GTP molecules. In this solution several processes take place 
(as well as their reverse reactions): 

(i) GTP hydrolysis: 

GTP — >GDP + Ai. (12) 

(ii) generating tubulin GDP from tubulin GTP: 

Tgtp — ► T G dp + A 2 , (13) 

(hi) growth of a MT: 

A 3 + MTjv-i + Tqtp — > MT N , (14) 

(iv) shrinkage of a MT: 

MT N — ► MT J v_ 1 + Tgdp + A 4 . (15) 

Note that experimental studies determined the values of the free energies for these reactions as: Ai ~ 
220 meV, A 2 — 160 meV and A3 ~ A4 ~ 40 meV, respectively [38]. These free energies are clearly 
above the thermal energy at room temperature (k&T ~ 26 meV) and they are within a quantum 
mechanical energy range that corresponds to the creation of one or a few chemical bonds. Hence we 
may consider each chemical reaction as a quantum mechanical process [39]. As a result, an individual 
microtubule with length L can be viewed as consisting of iV tubulin layers defining its quantum state 
\N). A tubulin layer consists of at least one tubulin dimer and at most 13 tubulin dimers as observed 
in the MT's structure. The state can be raised/lowered by a creation/annihilation operator (i.e. 
polymerization/depolymerization process) to the \N + 1)/\N — 1) state. The corresponding MT is 
then longer /shorter by one tubulin layer compared to the original one. 

In this paper, in order to simplify the problem we combine the above processes into two funda- 
mental reactions: 

(i) growth of a MT by one dimer by adding of one tubulin layer in an endothermic process: 

A + MTjv-i + T G tp — ► MT N , (16) 
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(ii) shrinkage of a MT by one dimer due to the removal of one layer of Tgdp dimer in an exothermic 
process: 

MTV — > MT^-i + T GDP + A, (17) 

where A is the energy of the reaction. In order to derive a quantum mechanical description of 
mechanisms (i) and (ii), we first need to introduce quantum states of MT, tubulin and heat bath: 

• \N) is the state of a microtubule with N dimers (both GTP and GDP tubulins). 

• \N T ) is the state of a tubulin, T GTP or T GDP . 

• | TV) is the GTP hydrolysis energy state. 

Then, the relevant second quantization operators would be: 



a t= \N + 1)(N\, 


(18) 


a=\N-l){N\, 


(19) 


&t = \N T + 1)(N T \, 


(20) 


b= \N T -1)(N T \, 


(21) 


dt= \N + 1)(N\, 


(22) 


d=\N-l)(N\, 


(23) 



Here b/tf and djS are annihilation/creation operators of tubulin and energy quanta, respectively. 
The operators a /a) are lowering/raising the number of tubulin layers that constructed a MT. Fol- 
lowing [40], one can express the above processes using creation and annihilation operators (18)-(23): 

a% d : A + MTV-i + T GTP — ► MI N (24) 
dH ] a : MTjv — ► MT w _i + T GDP + A (25) 

Operators (24) and (25) describe a MT's growth and shrinkage by one layer, respectively. Realis- 
tically, the polymerization or depolymerization process may happen repeatedly before reversing the 
process. This can be extended within our model by constructing product operators, i.e. (aJb d) m and 
(S a) n , where m and n are the number of growing or shrinking events in a sequence, respectively. 
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B. The Hamiltonian 



Based on the mechanisms in (24) and (25), the Hamiltonian for interacting microtubules with 
Tgtp/Tqdp tubulins can be written as 

k m l 

+ E M^k,- alb m d k _ m + A* m 4_ m a k ), (26) 

k,m 

where a;, -07, A and A are constants in units of energy. However, an intermediate transition between 
a microtubule in a growing phase and a microtubule in a shrinking phase must also be taken into 
account. A growing/shrinking microtubule may change its state quickly or after several steps to a 
depolymerizing/polymerizing state and then may change back to polymerizing/ depolymerizing state. 
Experimentally, the transition of microtubules from the growing to the shrinking phase is quantified 
by the catastrophe rate / cat and the transition from the shrinking to the growing phase is expressed by 
the rescue rate / res in which / rcs < / cat . As we discussed earlier, these transitions can be represented 
by a combination of creation and annihilation operators as the n th power of the reaction in (24) and 
(25): 



m 

oo 



+ E E ^[ A k n mJ n C k n mJ„ + A k n m n l n C l n m n I„]' ( 27 ) 

71=1 U ™ 1 

x^n 1 11 l n 1 l n — 1 



where 



c k n m n = (4i 6 mi d h )(al 2 b m2 d h ) . . . (ai n b mn d ln ). (28) 

Here k n = {ki, k 2 , . . . , k n } is a collection of indices and J] kn = Ski Sk 2 • • • Sk n - We note that the 
momentum conservation for the last two terms in the Hamiltonian (27) requires that 

n n n—1 

l n = ]Tk,-Em,-£L. (29) 

i=l i=l i=l 

Therefore, the first n — 1 of 1 will be free and summed in the Hamiltonian (27). 
In Bose-Einstein statistics the creation and annihilation operators satisfy 

[?k, qU = ^km, and [ql, q^} = = [g k , q m ], (30) 

where [A, B] = AB — BA is the Dirac commutator and q — a,b, and d . Since these operators 
mutually commute, the Eq. (28), can be rewritten as 

c k„iii„i n = a k!«k 2 • • • al n b mi b m2 ■■■b mn d h d h . . . d ln = a\b^ n d Xn , (31) 
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where l n is given by Eq. (29). 



IV. CLASSICAL EQUATIONS OF MOTION 

A system of coupled equations that describes the quantum dynamics of a MT is derived from the 
Heisenberg equation in Appendix A. However, MTs are overall classical objects (although some of 
their degrees of freedom may behave as quantum observable). As a result, we need to ensemble 
average over all possible states to obtain effective dynamical equations. 

Fourier transforming of a Vl b n and d v operators over all states, one can find 

V>(r, t) = n- 1 ' 2 ]T exp(-wj • r)a v (t), (32) 
v 

X (r, t) = 1T 1/2 J2 expHrj • r)b v (t), (33) 
v 

0(r,t) =fi- 1 / 2 ^exp(-^-r)rf T? (t), (34) 
v 

where Q is the volume over which the members of the plane wave basis are normalized [27]. Here 
■0(r, t), x( r , t), and 0(r, t) are the corresponding field operators for the quantum operators a v , b v and 
d v , respectively. The derivation of the equation of motion for the field operators is given in Appendix 
B. The final form of the equations of motion is found to be 

1 oo 

idf^ = A^ + iAx • W - ^^ 2 V 2 ^ + J2 (4™V*) r n ~\ n <p n , (35) 

^ n=2 
1 oo 

id tX = B oX + ■ V X - -B 2 V 2 X + E (S^V)^"- V B_ y (36) 

1 n=l 

i oo 

id t <P = C <P + id ■ - -C 2 V 2 + £(Cftyty B -V> +B "\ (37) 

1 n=l 

where n represents the degree of nonlinearity corresponding to the order of chemical process under- 
lying the term in the equation and Bi and C{ are constants which are given in Appendix. We 
obtain the general equations of motion for the system in terms of coupled nonlinear partial differ- 
ential equations (PDE's) that describe the MT field, the tubulin field and GTP field, respectively. 
These fields are complex function of space and time and their modulus squared corresponds to the 
spatio-temporal concentration of each of the three chemical species. 

In this paper we are primarily interested in the dynamics of MTs. Following Eq. (35), the dynamical 
equations for growing and shrinking states of a MT up to n = 3 can be written as 

id t i> = A Q ^ + iA 1 ■ - ^A 2 V 2 ^ + (4V) X 2 2 + (4V) <A W- (38) 
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Furthermore, the dynamics of the tubulin field, x, an d energy of the system, 0, are also determined 
by 

idtx = B oX + iB, ■ V X - \b 2 V 2 x + (39) 
id t 4> = C o + id • V0 - l -c 2 v 2 (t> + (cf (40) 

where, for simplicity, we just keep the n = 1 term. 

Although the system of equations (38)- (40) is generally similar to the phenomenological system 
of equations (3)- (6) found by other investigators, there are some manifest differences between them. 
Equations (38)-(40) are truly 3-dimensional that predict the most possible (or stable) 3D structure 
for a MT which is a vortex filament (see Sec. V) as observed experimentally. Other observed 3D 
structures of MTs such as double-wall MT, ring-shaped, sheet-like, C-shaped and S-shaped ribbons, 
and hoop structures (see [41-43] for details) during tubulin nucleation can be also explained through 
the model described here. Furthermore, the complex nonlinearity of MTs' behavior is embedded in 
our system through the role of fluctuations such as thermal noise that can generate catastrophe (see 
Sect. VB). 

As an example, the catastrophe event that is usually inserted by hand in those phenomenological 
models is a direct result of MT's fluctuations in our model (see Eq. 50). 

A vast array of mathematical methods of finding solutions to Eqs. (38)- (40) can be found in 
the monograph by Dixon et al. [44]. Among the analytical solutions known one can expect to find 
localized (solitonic) and extended (traveling wave) solutions. The latter ones may have the meaning 
of coherent oscillations that have been observed experimentally for high tubulin concentrations by 
Mandelkow et al. [45]. The localized solutions may correspond to nucleation from a seed. Sept et 
al. [46] studied the kinetics of a set of chemical reactions occurring during MT polymerization and 
depolymerization processes and also found oscillatory solutions with damping in a longer time regime 
of the assembled tubulins. 

V. THE DYNAMICAL EQUATION 

In order for the MT's assembly process to go through, one would expect that the energy in the 
system should be distributed uniformly during the course of experiment. This is based on the fact 
that one would expect the overall energy during the polymerization/depolymerization process to 
remain conserved. This, in fact, has been proved experimentally by direct measurements that the 
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enthalpy change during tubulin polymerization is about Oil kcal/mol of tubulin dimer [47]. It 
is interesting to note that in some experiments, the GTP-tubulin concentration in the solution is 
maintained to be constant, i.e. x = and V% = 0. The uniform energy distribution requires that 
the energy distribution function satisfies = and V0 = 0. As a result of this assumption, Eq. 
(40) can be solved for as 

~ -q i) X \ (41) 

where q is a constant parameter. Here and thereafter, we ignore the spatial derivative term in 
coefficients A iy Bi and Cj to avoid further complexities. Inserting Eq. (41) into Eqs. (38) and (39) 
we find 

id t ip + iv- Vip = --bV 2 ^ + Vifj, (42) 

td tX = fV 2 x + U X , (43) 

\/(|^|,|x|) = a + c|x| 4 |^r-rf| X | 6 |^| 4 , 
f/(H)=e-/#| 2 , 

where \ip\ 2 = ipip* . In the above equations we introduced a set of new parameters for simplicity. We 
also transformed r to r + and then set Ai — = v where v can be considered to represent a 
MT's velocity relative to the tubulin concentration in the solution. Here, parameters a, b, e and / are 
real but c, d and h are complex. Eq. (42) represents the nonlinear cubic-quintic Schrodinger (NLS) 
equation with a complex potential that has been extensively studied in connection with topics such as 
pattern formation, nonlinear optics, Bose-Einstein condensation, superfluidity and superconductivity, 
etc. [48]. In a series of papers, Gagnon and Winternitz discussed symmetry groups of the NLS 
equation and provided some exact solutions in spherical and cylindrical coordinates [49-53] which is 
of relevance to the present case. General solution of the NLS equation can be cast in the form of 
■0(r, t) = R(r, t) exp[iS(r, t)] which involves topological defects (point in 2D and line in 3D). In three 
dimensions these defects represent one-dimensional strings or vortex filaments [48] . In a cylindrical 
coordinate system, there exists a stationary solution that represents a straight vortex filament with 
twist: 

ip(r, 9, z, t) = R(r) exp{i(ut + n6 + w(r) + k z z)}, (44) 

where u> is the spiral frequency, R(r) is the amplitude, w(r) is the spiral phase function and integer n 
is the winding number of the vortex [54]. The axial wave number k z characterizes the vortex's twist. 
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k z = represents an untwisted vortex that is the most stable solution [48]. In the case of the NLS 
equation, a family of vortices that move with a constant velocity is also a solution [54]. 



A. Numerical results 

Equations (42) and (43) are solved numerically with a no-flux boundary condition. As an initial 
condition we chose a straight vortex filament perturbed by small noise (eg. thermal or environmental 
noise). In Fig. 1 we compare the observed data on the MT length as a function of time with our 
simulation results. The length of a vortex is defined as [54] : 

L(t) = je^o-hAM)|)tfV, (45) 

where Q(x) is the step function and ipo is a constant. In Figs. 1 and 2 we compare the observed 
data on the MT length as a function of time with our simulation results. Experimental panels in 
Figs. 1 and 2 represent the experimental data published by Rezania et al. [55]. Simulation panels 
show the numerical results of the normalized vortex length as a function of time for the given set of 
parameters. 

To provide a simple yet accurate and powerful comparison between experimental and simulation 
results, we graph recursive maps for the data points. The advantage of the recursive maps is the 
introduction of regularity into the data sets that allows for a better choice of adjustable parameters 
due to noise reduction inherent in the separation of data into subsets corresponding to independent 
processes. In spite of being very simple, recursive maps of assembly and disassembly processes of 
individual MTs can successfully reproduce many of the key characteristic features. Consider first the 
following stochastic map as the simplest case that illustrates the approach taken: 

£(t n + l)=r[£(t n ) + 5], (46) 

where £(t n ) is the length of a microtubule after n time steps, t n . The parameter r is chosen to be a 
random number with the following two possibilities: 




1 with probability p 
with probability 1 — p 



In terms of the MT polymerization process, p is the probability that a given event will result in 
assembly while 1 — p is the probability of a complete catastrophe of the MT structure. The above 
simplified model, therefore, is governed by only two adjustable parameters: (i) the probability of a 
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complete catastrophe 1 — p which is constant and independent of the length or time elapsed and 
(ii) the rate of polymerization which is proportional to the length increment 5 over the unit of time 
chosen in the simulation. Thus, the coefficient 5 divided by the time step At (= t n+i — t n ) gives 
the average growth velocity of an individual MT. Such information can be used to fine tune the 
simulation parameters. We note that the slope of the line in the simulation panels in Fig. 1 can 
be adjusted by varying the real parts of parameters c and d. The frequency of catastrophe events 
can also be changed by adjusting the parameter b. In the recursive map panels in Figs. 1 and 2 we 
compare the recursive map for both the experimental data and the simulation results. Based on the 
recursive maps, the key characteristics of the experimental and simulated results that were obtained 
independently are quite similar. This represents Eqs. (42) and truly describes the dynamics of MTs' 
polymerization. 

To provide a more solid comparison, a spectral analysis is also carried on both experimental and 
simulated data. As discussed by Odde et al. [56], the power spectrum analysis is a more general 
way to characterize the microtubule assembly /disassembly dynamics without assuming any model a 
priori. The power spectrum panels in Figs. 1 and 2 represent the spectral power of the experimental 
and simulated data, respectively. As shown, there is a great agreement between the experimental 
(solid curve) and simulated (dashed curve) spectrums. We note that the curves are plotted at different 
offsets for visual clarity only. As expected, no particular frequency of oscillation can be found from 
the power spectrums. However, both power spectrums demonstrate a very similar broad distribution 
that more or less decays with frequency as an inverse power-law with slope ~ 1.0 in Fig. 1 and ~ 1.2 
in Fig. 2, respectively. The best fit inverse power-law is shown by a dotted line in both panels. 

More interestingly, our results show no attenuation states during MT polymerization (Fig. 1). 
The MT length undergoes small fluctuations all the time. This can be understood by noting that 
our model is based on the cyclic polymerization and depolymerization of tubulin dimers. Behavior 
consistent with this result has recently been observed by Schek et al. [57] who studied the microtubule 
assembly dynamics at higher spatial (~ 1-5 nm) and temporal (~ 5 kHz) resolutions. They found 
that even in the growth phase, a MT undergoes shortening excursions at the nanometer scale. 
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B. Analytical results 



Using Eq. (45), the time- varying nature of vortex filaments' length can also be extracted analyt- 
ically. Taking time derivative of L(t) one finds 

= ~/r 1 i 

2 J 

= ~\ j j^j [-v • V W) + (i6/2)(^*VV - ^V 2 f) - i(V - O^Wo - IVDA 

where £ = ^ — IVK 1 ", t)\ and is the Dirac delta function. For a general vortex solution, this may 
lead to a very complicated function of time. However, for the given vortex solution, Eq. (44), the 
above statement will be simplified as 

^-L(t) ~ f[v ■ VR + (6/4) (4 VR ■ Vw + 2RV 2 w) + Im(V)R}5(tp - \tp\)d 3 r, (47) 
at J 

where the above integral is constant in time. As a result, the vortex length can grow in time linearly. 
In order to calculate possible fluctuations in the vortex length, we introduce a small perturbation 
(due the thermal or environmental noise) 5ip(r,t) ~ £(r) exp(icu't) where |£| <C 1 and uj' represents 
the quasi-periodicity of the fluctuations. Inserting it into Eq. (45) and taking the time derivative, 
we have 



d_ 
~dt 



L + 5L)= J d t Q{^o - |V + S^\)d 3 r, 



I j U + tw/ ^* + + WWlfo - \^\)d 3 r. (48) 



2 J |^ + #| 

Taking the time derivative and expanding the denominator in the integral we have 

dtiw*) = -v • V(W*) + it[V*V 2 V - ^vV - i(V - v*)w% 

d t (ip6ip* + tp*6ip) = -v • Vip Sip* + -Sij*V 2 ij - i{u' + V)ip5ip* + c.c. 

2 



m 



1 - lw* sr - lw rs * 



where c.c. means the complex conjugate. Again, for the given solution (44) one finds after some 
manipulations 

r 

5L^i— sin(cj — uj')t + const., (49) 

D — UJ 
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where 

T = y*[- v • VF + (*/2)V 2 F - i(y + u)')F - (l/2)F/|F| 3 ]r 5(^ - M) d 3 r + c.c, (49) 

and F(r) = R(r) exp[i(w(r) +116 + k z z)\. Interestingly, from Eq. (49) one can see that even if u/ 
goes to zero, the vortex length fluctuates with the spiral frequency uj. As a result, the vortex length 
fluctuates all the time. Furthermore, when lu' — > uo Eq. (49) reduces to 

SL fa -Tt + const., (50) 

where represents a linear shortening similar to catastrophe events. As a result, a catastrophe event 
can be explained within our model. 

VI. DISCUSSION 

In our model, the basic structural unit is the tubulin dimer. Each dimer exists in a quantum 
mechanical state characterized by several variables even in our simplified approach. Each microstate 
of a tubulin dimer is sensitive to the states of its neighbors. Tubulin dimers have both discrete 
degrees of freedom (distribution of charge) and continuous degrees of freedom (orientation). A 
model that focuses on the discrete ones will be an array of coupled binary switches [58, 59], while a 
model that focuses on the continuous ones will probably be an array of coupled oscillators [7, 60]. 
In the present paper we have focused on tubulin binding and GTP hydrolysis as the key processes 
determining the states of microtubules. These are also the degrees of freedom that are most easily 
accessible to experimental determination. In this paper we have shown how a quantum mechanical 
description of the energy binding reactions taking place during MT polymerization can lead to 
nonlinear field dynamics with very rich behavior that includes both localized energy transfer and 
oscillatory solutions. 

In particular, based on the chemical binding reactions that are taking place during microtubule 
polymerization, a quantum mechanical Hamiltonian for the system is proposed. Equations of motion 
are then derived and transformed from the purely quantum mechanical description to a semi-classical 
picture using the inverse Fourier transformation. After lengthy calculations we found that the dy- 
namics of a MT can be explained by the cubic-quintic nonlinear Schrodinger equation (NLS) with 
variable coefficients. A generic solution of the NLS equation in cylindrical geometry is a vortex fila- 
ment [48, 53, 54]. Interestingly, we showed both analytically and numerically that such a solution can 
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grow or shrink linearly in time as well as fluctuate temporally with some frequency. This behavior 
exhibits two distinct dynamical phases: (a) linear growth/shrinkage and (b) oscillation about some 
mean value, and is consistent with the characteristics of the MT's dynamics as observed in different 
controlled experiments in vitro (Fig. 1). 

It is noteworthy that dynamics of pattern formation can be also described by NLS equation in 
which ip(r, t) represents the order parameter. Interestingly, a number of convincing experiments, 
performed by Tabony et al. [61] demonstrated that gravity can indeed influence certain chemical 
reactions. Tabony and his colleagues, at the French Atomic Energy Commission lab in Grenoble, 
found that when cold solutions of purified tubulin and the energy-releasing compound GTP were 
warmed to body temperature, microtubules formed in distinct bands. These bands form at right 
angles to the orientation of the gravity field or, if spun, to the centrifugal force. Despite several 
studies [62-65], the above experiments are yet to be fully explained theoretically. Our goal in future 
studies is to focus on the dynamics of pattern formation by MTs using the results presented in this 
paper. 

We have demonstrated here that the assembly process can be described using quantum mechanical 
principles applied to biochemical reactions. This can be subsequently transformed into a highly 
nonlinear semi-classical dynamics problem. The gross features of MT dynamics satisfy classical 
field equations in a coarse-grained picture. Individual chemical reactions involving the constituent 
molecules still retain their quantum character. The Fourier transformation allows for a simultaneous 
classical representation of the field variables and a quantum approach to their fluctuations. Here, 
the overall MT structure (and their ensembles) can be viewed as a virtual classical object in (3+1) 
dimensional space-time. However, at the fundamental level of its constituent biomolecules, it is 
quantized as are true chemical reactions involving its assembly or disassembly. 
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APPENDIX A: DERIVATION OF THE HEISENBERG EQUATIONS OF MOTION 

The Heisenberg equation of motion for a space- and time- dependent operator q(r, t) reads as 

ihdtq(T,t) = -[H,q(T,t)], (Al) 

where H is the Hamiltonian. Before finding equations of motion, one needs to calculate the commu- 
tation relation [q„, qj- ] that is 

[?ti, <ij = ?k! • • ■ ?tj = ^,k! qi 2 qls ■ ■ ■ iL + s vm qIA 3 • • • i + • • • 

+^,k n qLA, ■■■<lL-l- ( A2 ) 

Since all ki, k 2 , . . . , k„ are dummy indices one can write Eq. (A2) as 

[<lr» Qt } = n ^,k„ Qt , , (A3) 

where k„ is chosen for simplicity. Using Eq. (A3) we can find the commutation relations between 
a n and b n operators with a- - ? and c~ ~ T operators as 

K C Lm„lJ = K 4„ & k«kJ = « ^m^^k.^ic^ (A5) 

However, the commutation relation between oL operator and c~ . T will be 

k , , IT1 n 1 , ( 

\d„,cl ~ T 1 = [d„, d~ at 1 = — 1) 5ni ,d~ d\ + 5„i dl ) b'~ at , (A6) 

where 1„ is given by Eq. (29). Therefore, the equation of motion for a v , b v and d v operators (and 
their Hermitian conjugates) can be derived from Hamiltonian (27) as 



idtdn = to n a n + n A^r - t a~ 6™ dr d ,v^n-i n , v^n , (A7) 

n k n _irii„l n _i 

i<9+&„ = vj „b „ + n Ar - ? d~ d™ v^n-i, *&Jh , (A8) 

t T7 ?7 T7 ^ T7k n m n _il n _i i n _ 1 ^)" =1 ki-»j-^)" =1 (mj+lj) m "- 1 k ™ ' V 7 

™ k„m„_il n _i 

id t d v = + E _ E (n - 1) Vnm^-^L.^^k.-mO-rj-En 2 U 

™ k„m„l„_ 2 

+ ^ ^ ^ Y^ n ci ^ v^n-l i Ar - T d| b~ Ofl . (A9) 



^nlUn-ln— 1 



19 



APPENDIX B: DERIVATION OF EQUATION OF MOTION FOR THE FIELD OPER- 
ATORS 



Multiplying both sides of Eq. (A7) by exp(— 277 ■ r), dividing by Vl 1 / 2 and summing over 77, one finds 



id t ijM- 112 



^c^expt-irj • r)a v 



+ y^y^ ^ nAc - t exp(— -m-r) a\ d~, d , v^n-i,. . . 

Z-^i Z-^i r7k n _im n l„_i r\ I J k„_i mn l n-i rj+ > (kj -lj) — V _ , : 

n *7 k rh 1 1 



(Bl) 



Changing 77 — > ry — S"=i (kj — 1?) + Z)"=i m j m the second term of Eq. (Bl), one finds 



id t ij = n- 1/2 



5^w„exp(-M7 • r)a„ 
v 

+ EE_ E n Atj-£ kn-imj^-i' 
n *? k n m n l„_i 



177-r+i J]™ j 1 (kj -lj )-r-i X;" =1 mj t 



(B2) 



or 



id t ijM~ 1/2 



E^ ex pH T 7- r N + EE E ^V^^^u 

K. n IIl n l n — 1 



T e -™,r & e -il«-iT rf 
k n _l m ™ ln-1 



(B3) 



where £ = X^=i(kj — 1?) — Z)^=i m i- Here, for example, exp(— ik n ■ r) = exp(— iki • r) exp(— ik2 • 
r) . . . exp(— ik n ■ r) = exp(— i£)" =1 kj ' r )- Our goal is now to rewrite Eq. (B3) in terms of field 
operators, tp, Xi 4>i an d their derivatives. This can be done in a straightforward manner provided 
the dispersion matrix elements u v and A ^ k„_im„I n _i which are generally function 77, kj, nij and 
h (1 < i < n) are known. Unfortunately, such information is very model dependent. Therefore, the 
simplest way that also keeps the generality of the problem is to Taylor expand these matrix elements 
about some point (ry , k j, m j, loi) in the space spanned by 77, kj, rrij and lj [27, 36]. 
Expanding u v to all orders, one finds 

00 

Wr, = ^0 + E " »7o) ' V^l'wo/s! , (B4) 

s=l 

where u = uj Vq . Furthermore, for any function f(rj, k„, m„, l n ) = A j^^j^ we can write 

/(77,k n ,m n ,l n ) =f + (n- rio) ■ V v f\o + E( k i ~ k oj) ' V k ./| 

3=1 
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E( m , - m ,) • V m ,/| + E& " ■ Vij./lo + 

3=1 3=1 



n oo s u s—u—v 



+ E EEE E sc » T » m c w / s \ 

p,q,r=l s=2 u=0 v=0 w=0 

x [fa - T7 ) • V,,]" [(k p - k 0p ) • V kp f [K - m 0g ) • V m J^ [(l r - l 0r ) ■ V lr ] s - u - v - w f\o, 

(B5) 

where 

n 

fo= E /(^o^op^o^W), (B6) 

p,q,r=l 

where s C r are binomial coefficients. Here, for example, V m / means id mx f + jd my f + kd mz f where 
i,j and are unit vectors in the m x ,m y and m z directions, respectively, and V m /| is the value of 
the gradient at point (r} , k , m , lo). 

Using Eqs. (B4) and (B5), Eq. (B3) can be written as 

id t i> = a (u# + iMu) ■ vv - \ YXM^hdl^ + E r^ 3 "- 1 )/ 2 AS n V +n_ y<r 

i,j n 

+ J2nn^-^ 2 U^-'x^-'v^io • V0 + E v k ,/| • w + i> +n ~\ n r 

n \ j=l 

n n—1 \ 

+ E^ + " _1 V m ,/| • V X X-V" + E ^""VVi./lo • V0 0- 1 , (B7) 

J'=l J'=l / 

where 

Ao(w) = cj - t/ • V„w| + (1/2) E^ojdJ^lo , (B8) 
[Ai(w)]i = - E^-^.^lo + <9^| , (B9) 
[A 2 (o;)]« = ^"lo, (BIO) 

n— 1 n 1 

A?° = /o - »?o ' V,/|o - E ko; • V fe ,/|o - E m o, • V m3 /| - E hi ■ V^/lo. (Bll) 

i=i i=i j'=i 

Similarly, using Eqs. (A8) and (A9), one can write equations of motion for x an d as 



i^X = AoMx + «Ai(c7) • V X - \ YXM^l^X + £nfi< te - 1 >/ 2 A<Vx +n -V +n 

+5>n< 3 »- i >/ 2 (Vx+^V+^v^/io • V0+ + e v k ,/| • w r~ l x +n <P +n 

n \ j=l 

n—1 n—1 \ 

+ E ^ n V mj /|o • V X + X n " 2 0" + E ^V"" 1 V^/lo • V0+ , (B12) 

J'=2 J=l / 
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id^ = X (a)cf> + i\ 1 (a) • - \ E[A 2 (a)]^0 + £ ^"-^A^V"^"" 1 



n 



+ £ ^ (3n - 1)/2 (n - l)^x + > +n_2 V^|o • V0 + + n E V kj /| • ^""V V" 1 



3=1 

n— 1 



3=1 3=1 



(B13) 



where 

n n— 1 n— 1 

A 2 n) = /o - i? • V„/| - E k o, • V fcj /| - E m o, • V m ./| - E V- • Vj^/lo, (B14) 

i=i i=i j'=i 

7i ?2 n. — 1 

Aj° = /o - rjo ■ V„/| - E k o, • V fcj /| - E m o, • V m , /|„ - E W • V^/| . (B15) 

J=l 3=1 3=1 

Simplifying the equations of motion as 

id t i> = A *p + iA, ■ - \a 2 V 2 ^ + E (4V) ^ +n ~\ n r, (B16) 
id tX = Box + iBi ■ V X - \b 2 V\ + E(^ n) W^V"" 1 ^, (B17) 



2 

id t <f> = C o + id • V0 - ic 2 V 2 + E(^ n V)^ n - 1 X + > + "" 1 , (B18) 



where 

n-1 

A = AoH, Ai = Ai(w), A 2 = \ 2 (lu), 4 n V + = n^(Ai n) + EV kj /|o-V)^+ 

(B19) 

B = Xo{w), Bi = Ai(ro), 5 2 = A 2 (ro), ^ n V = nfi^(A^ + EV kj /| - Vty, 

i=i 

(B20) 

n 

C = \ (a), C^X^a), C 2 = \ 2 (a), Cf V = ^^(A^ + E V k ,/| • V)V>. (B21) 
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TABLE I: The parameters available in the literature. 



Parameter Simulation Coeff. Exp. Value Reference 

MT growth rate Real(c) 0.50 - 19.7 (/xm/min) [20, 66-72] 

MT shortening rate Real(d) 4.1 - 34.9 (/xm/min) [66-68, 70-72] 



MT catastrophe frequency 0.12 - 3.636 (/min) [67, 68, 70, 71] 

MT diffusion constant b 2.6 - 30.3 (^m 2 /min) [73-75] 

Tubulin diffusion constant e 300 - 480 (>m 2 /min) [12, 76] 
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Figure Legends 



Figure 1.: Length of a distinct microtubule as a function of time. The top-left panel represent 
experimental data published in [55]. The top-right panel is the simulation result with the set of 
parameters (a = 1, b — 10, c = 10 + i, d — 20 + i, e = 300, / = 1 and h = —.1 + i). The bottom- 
left panel shows recursive maps for both experimental and simulation results. The bottom-right 
panel represents the power spectrum of the experimental (solid curve) and simulated (dashed curve) 
data, respectively. The curves are plotted at different offsets for clarity. No particular frequency 
of oscillation can be seen from the power spectrums. However, both power spectrums show a very 
similar broad distribution that more or less decays with frequency as an inverse power-law with slope 
~ 1.0. The best fit inverse power-law is shown by a dotted line. 

Figure 2.: Same as Fig I. but with the set of parameters (a = 1, b = 30, c = 10 + 10i,d = 
20 + lOi, e = 300, f = 1 and h — — .1 + i). In the power spectrum panel, the inverse power- law has a 
slope of ~ 1.2. 
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